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Abstract 


A  set  of  second  order  elliptic  partial  differential  equations 
for  the  generation  of  three-dimensional  curvilinear  coordinates 
between  two  arbitrary  shaped  bodies  have  been  proposed*  The  resulting 
equations  have  only  two  independent  variables  and  therefore  require 
an  order  of  magnitude  less  working  core  capacity  than  when  the  equa¬ 
tions  depending  on  all  the  three  independent  variables  are  considered. 

The  resulting  equations  have  been  programmed  for  the  numerical 
solution  of  the  equations  on  the  Cray-computer .  Much  of  the  time  has 
been  spent  on  the  generation  of  surface  coordinates  for  an  aircraft 
fuselage  by  spline  and  various  other  methods.  Once  these  surface 
coordinates  have  been  accurately  established,  the  proposed  field  equa¬ 
tions  will  be  solved  for  the  region  between  the  fuselage  or  other  body 
and  another  arbitrarily  selected  outer  surface  (e.g.,  a  sphere).  The 
spline  method  of  the  body-coordinates  is  discussed  below. 

It  has  been  established  that  spline  interpolation  can  be  used 
to  construct  a  computational  grid  about  a  simple  wing-fuselage  con¬ 
figuration.  For  generality,  the  components  of  the  configuration  are 
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A  Method  for  the  Generation  of  General  Three-Dimensional 
Coordinates  Between  Bodies  of  Arbitrary  Shapes* 

by 

Z.  U.  A.  Warsi^ 

Department  of  Aerospace  Engineering 
Mississippi  State  University 
Mississippi  State,  MS  39762 

Abstract 

Analytical  development  of  a  set  of  second  order  elliptic  partial 
differential  equations  for  the  generation  of  three-dimensional  curvi¬ 
linear  coordinates  between  two  arbitrary  shaped  bodies  is  presented. 

The  resulting  equations  have  only  two  independent  variables  and  therefore 
require  an  order  of  magnitude  less  working  core  capacity  than  when  equations 
depending  on  all  three  independent  variables  are  considered.  The  method 
also  allows,  in  a  straight  forward  manner,  the  possibility  of  coordinate 
contraction  in  the  desired  regions. 

An  exact  solution  of  the  proposed  equations  for  the  case  of  an  inner 
prolate  ellipsoid  and  an  outer  sphere  with  coordinate  contraction  is  pre¬ 
sented  to  demonstrate  that  by  using  these  equations  it  is  possible  to 
generate  three-dimensional  coordinates  between  analytically  specified 
surfaces  of  simple  forms  by  analytical  means. 

The  fundamental  constraining  equations  which  have  been  adopted  for 
the  generation  of  coordinates  are  =  0  and  A2h  =  0,  where  is  the 
surface  Beltrami  operator  of  the  second  order. 

♦Research  supported  in  part  by  the  Air  Force  Office  of  Scientific 
Research,  under  Grant  AFOSR  No.  80-0185. 
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1.  Introduction 


At  present  a  number  of  techniques  are  under  active  development 
for  the  generation  of  three-dimensional  body-oriented  coordinate  systems 
for  use  in  the  numerical  solution  of  the  Navier-Stokes  equations  and 
other  field  equations  where  the  exact  specification  of  the  boundary 
conditions  is  of  prime  importance.  Among  these  efforts  two  easily  dis- 
cernable  groups  can  be  formed,  (i)  algebraic  methods,  and  (ii)  the 
elliptic  equations  method.  In  the  first  group  the  grid  points  in  space 
are  obtained  by  some  interpolation  or  blending  functions  scheme  which 
depends  on  the  given  boundary  data.  The  choice  of  the  interpolation 
scheme  or  of  the  blending  functions  is  crucial  in  achieving  a  desired 
order  of  smoothness  and  distribution  of  the  grid  points  in  space.  This 
line  of  effort  has  actively  been  considered  by  Eiseman  [1,2],  Smith  and 
Weigel  [3],  and  Eriksson  [4],  In  the  second  group  of  efforts,  a  set  of 
three  poisson  equations  in  the  curvilinear  coordinates  are  first  inverted 
and  then  solved  for  the  Cartesian  coordinates  under  the  prescribed  values 
at  the  given  boundaries.  Thus  in  essence  all  the  methods  of  the  second 
group  are  a  straight  forward  extension  of  the  work  of  Thompson  et  al  [5] 
in  two  dimensions.  Research  in  this  area  has  been  conducted  by  Mastin 
et  al  [6],  Yu  [7],  Ghia  et  al  [8],  and  Graves  [9]. 

At  this  stage  of  research  it  is  premature  to  compare  the  two  groups 
since  neither  of  them  have  been  fully  investigated  for  their  inherent 
potentials.  However,  based  on  the  success  of  the  differential  equations 
approach  in  two  dimensions,  e.g.  [5],  it  is  desirable  to  further  investi¬ 
gate  the  elliptic  equations  approach  for  the  generation  of  coordinates. 

The  elliptic  equations  approach  presented  in  this  paper  is  different 
from  the  approaches  adopted  in  the  previously  cited  works,  i.e..  References 
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[6]  -  [9J.  The  proposed  method  depends  heavily  on  the  formulae  of 
Gauss  and  on  the  concept  of  principal  curvatures  of  a  surface.  It  has 
been  shown  that  a  fruitful  arrangement  of  the  classical  differential- 
geometric  results  can  yield  a  method  which  is  easily  programmable  on  a 
computing  machine,  and  which  at  any  time  solves  a  two-dimensional 
partial  differential  equation  of  the  form  used  in  Ref.  [5].  In  this 
paper  only  the  theoretical  development  of  the  method  along  with  a  tech¬ 
nique  to  redistribute  the  coordinate  surfaces  near  the  inner  boundary 
surface  has  been  considered.  The  developed  equations  have  been  solved 
for  the  generation  of  three-dimensional  coordinates  between  an  inner 
prolate  ellipsoid  and  an  outer  spherical  surface  in  an  exact  analytic 
form. 
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2.  Notation  and  Collection  of  Formulas 


In  what  follows,  the  general  coordinates  are  denoted  as  x*  (i  =  1, 
2,3).  However  when  an  expression  has  been  expanded  out  in  full  and  there 
is  no  use  for  an  index  notation,  we  have  set 

x1  =  £  ,  x2  =  n  ,  x3  =  £  . 

The  derivatives  of  the  position  vector  r  =  (x,y,z)  are  denoted  as 


r . 
~  i 


*ij 


a2r  _ 

3x*3x^ 


The  covariant  components  of  the  metric  tensor  are 


8ij  = 


while  the  contravariant  components  are  given  by 


ij  _i 

g  8kj  =  6k 


Thus  in  three  dimensions 


(2.1) 


(2.2) 


g  =  det  (g-^j) 


8ll822g33  +  2gl28l3823  ~  ^823^2gll  ~  ^813^2g22  “  ^812^2g33 


(2.3) 


Writing 


G1  =  g22g33  ~  ^g23^ 2 
G2  "  gllg33  ~  (g13)2 
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G3  =  811822  “  ^g12^2 
G4  =  813823  “  812833 


G5  =  812823  "  813822 
G6  =  812813  "  gllg23 


we  have 


g11  =  Gl^g  »  g22  =  G2^ 8  »  g33  =  G3^ g 

g12  =  G4/g  ,  g13  =  G5/g  ,  g23  =  G6/g 


The  Christoffel  symbols  based  on  the  metric  g 


ij 


are 


rjk  =  gi£[jM) 


where 


[jk,£]  =  rjk 

■*k,  * 

2\  k  +  „  1 


i  > 
3x 


(2.4) 


(2.5) 


(2.6) 


and  repeated  lower  and  upper  indices  imply  summation.  In  the  sequel  we 
have  also  used  the  surface  Christoffel  symbols  which  have  been  denoted  as 

,  where  the  Greek  indices  range  over  (1,2)  or  (3,1)  or  (2,3). 

PY 

The  coordinate  which  is  held  fixed  to  account  for  the  surface  geometry 
is  denoted  by  a  superscript  in  parentheses.  Thus  the  unit  normal  vector 
on  the  surface  v  ■  const .  is  given  by 


5 


mmi 


(2.7) 
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V 


r 

~aa 


aS 


=  r 
aa~y 


=  TY  r 

a3~y 


=  TY  r 
SB-y 


+  S 


+  T 


+  U 


(v)n(v) 

(v)n(v) 

(v)n(v) 


(2.12) 


Where  (v,a,8)  are  in  the  permutational  sequences  of  (1,2,3)  as  shown  in 
(2.8),  and  the  repeated  index  y  implies  summation  on  the  two  indices  of 
a  surface. 

The  sum  of  the  principal  curvatures  of  the  surface  v  =  const,  is,  [10], 


t(v)  4.  i<v)  -  (n  n<v)  o„  4-  0  e<v)Wr 

kl  +  k2  "  (8aaU  "  2gaBT  +  gBBS  )/Gv 


(2.13) 


where  in  writing  equation  (2.13)  for  a  particular  value  of  v,  use  must  be 
made  of  Eqs.  (2.8)  and  (2.10).  We  now  introduce  two  second  order  surface 
differential  operators  by  using  (2.8),  which  for  v  =  const,  are 


D  "  8B38aa  ~  2gaB3aB  +  gcta8gB 


(2.14) 


A2V)  z  |^-l3a{^g-(8eeao  ~  ga$V} 

v  v 


+  V^aa’l!  -  8«eV» 

V 


(2.15) 


As  is  well  known,  the  operator  the  Beltrami  differential  operator 

of  the  second  order  [10]. 

The  three  space  Christoff el  symbols  which  have  been  referred  to  in 


the  next  section  are  given  below. 


3.  Formulation  of  tho  Problem 

The  principal  idea  of  the  method  to  be  presented  is  to  generate  a 
series  of  surfaces  on  each  of  which  a  certain  a'priori  chosen  variable 
or  coordinate  is  kept  fixed.  Each  surface  to  be  generated  starts  from 
a  given  curve  of  the  inner  body  and  ends  on  the  corresponding  curve  of 
the  outer  boundary,  cf.  Fig.  1.  A  routine,  preferably  a  spline  fit, 
can  then  be  used  to  join  the  successive  generated  surfaces  so  as  to  have 
a  smooth  three-dimensional  computational  net  for  solving  other  physical 
field  equations. 

To  illustrate  the  method,  we  take  v  =  3  or  x3  =  £  =  constant  on  each 
surface  to  be  generated.  Thus  a  =  1  and  8  =  2,  viz.,  a  and  8  respectively 
correspond  to  the  coordinates  x1  =  £  and  x2  =  n-  For  the  sake  of  brevity 
of  notation  we  will  not  use  the  superscript  (3)  unless  it  becomes  neces¬ 
sary.  Thus  from  (2.10) 

S  =  n  •  r,_  ,  T  =  n  •  r_  ,  U  =  n  •  r  (3.1) 

~£n  -  -nn 


where 


n  =  iX  +  jY  +  kZ 


(3.2) 


Equations  (2.12)  are 


hi  '  tIiey  +  85 


(3.3) 


r  =  T'  r  +  Tn 
-£n  12-y 


(3.4) 


r 

-nn 


T,,r  +  Un 
22~y 


(3.5) 
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(31  (31 

From  (2.14)  and  (2.15)  the  operators  Dv  and  4^ 


are 


D  H  8oo  — o  “  28 


a2 


22  ..2  “e12  a^an  +  8n  .  2 


(3.6) 


U 


an 


<  -  l  f a  f  i  /  a  3  ii 

a2  =  >^-[ae{>^-(s 22  as  8i2  an)} 


,  a  f  l  ,  a  a  v  ■>  •. 

+  an^^g-^n  an  “  812  ac 


(3.7) 


We  now  multiply  equations  (3.3)  -  (3.5)  respectively  by  ~28i2’ 

g  adding  and  using  equations  (2.13)  and  (2.15)  to  have 


Dr  +  G3(r^A2£  +  r^A2n)  =  G^n(k1  +  k2) 


(3.8) 


where 


A2^  G3(2g12T12  ~  822T11  811T22^ 


A2n  G3(2g12T12  822T11  "  811T22) 


(3.9) 


G3  ”  811822  ~  (s12^2 


To  obtain  an  expression  for  k^  +  k2  consider  equation  (2.11)  and 
utilize  the  property  that  n  is  orthogonal  to  r^  and  r^,  so  that 


5  '  '  rh<5  ‘  5t> 


ru<!>  •  Ei> 
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n  •  r  =  rL(n  •  r  ) 

-  ~nn  it-  - 

where  all  the  derivatives  with  respect  to  £  are  evaluated  at  ?  = 
constant.  Multiplying  Eq.  (3.8)  scalarly  by  n  and  using  (3.10), 
get 


(3.10) 


we 


G3(kL  +  k2)  =  (n 


rc)(g 


-  2g 


11  22 


12ri2  +  822ril^ 


(3.11) 


We  now  propose  the  following  deterministic  problem:  Let  £  and  n  be 
the  surface  coordinates  on  the  surface  £  =  constant,  subject  to  the 
constraints 


A2£  =  0 
A2n  =  0 


(3.12) 


Then  the  Cartesian  coordinates  x, y,z  of  the  surface  satisfy  the  differential 
equations 


Dr  =  G3(k3  +  k2)n 


(3.13) 


The  three  scalar  differential  equations  for  the  generation  of  the  Cartesian 
coordinates  are  then 


*22*55  '  2812*5n  +  ellxnn  *  XR 
822”{5  -  2g12y5n  +  8llynn  ’  YR 
*22z«  -  2g12zCn  +  8llznn  *  ZR 


(3.14) 

(3.15) 

(3.16) 


where 
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(3.17) 


R  =  (Xx  +  Yy  +  Zz?)(811I’ 22  ~  2g12ri2  +  822F11^ 


and 


X  -  -  Vc)/’^3 


V  ‘  <V{  -  XE2r,)//5I 


(3.18) 


z  -  -  xnyc)//^ 


Equations  (3.14)  -  (3.16)  form  a  quasilinear  system  of  partial 
differential  equations  in  which  the  components  of  r^  are  assumed  to  be 
known.  Since  the  values  of  x,y,z  are  known  on  the  basic  inner  and  outer 
boundaries  (denoted  at  B  and  »  respectively  in  Fig.  1),  a  suitable  way 
of  prescribing  r^  can  be  to  take 


h  =  fl(n)(~c}B  +  f2(n)(^)o 


(3.19) 


where  f^(n)  and  f2(n)  are  suitable  weights  having  the  properties 


W  ‘ 1  ’  f2(V  - 0 


£l(n.)  '  0  '  £2(n«)  '  1 


For  exposing  the  essential  nonlinear  terms  in  the  factor  R  we  refer 
to  Eqs.  (2.16)  -  (2.18)  in  which  the  r^  terms  have  been  collected 
separately. 

Referring  to  Figure  2,  we  now  solve  Eqs.  (3.14)  -  (3.16)  for  each 
C  *»  const.,  on  a  rectangular  plane  by  prescribing  the  values  of  x,  y  and 
z  on  the  lower  side  (C^)  and  upper  side  (C^)  which  represent  the  curves 
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on  B  and  »  respectively.  The  sides  C3  and  are  the  cut  lines  on  which 
periodic  boundary  conditions  are  to  be  imposed.  The  preceding  analysis 
thus  completes  the  formulation  of  the  problem. 
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4.  Coordinate  Transformation  (Contraction) 


For  the  purpose  of  generating  coordinates  between  the  space  of  the 
inner  and  outer  boundary  which  can  be  distributed  in  a  desired  manner, 
we  consider  a  coordinate  transformation  from  £  X  and  n  cr.  Let 


e  =  ?(x)  +  ?0 


n  =  n(a)  +  nB 


(4.1) 


then 


Writing 


«<*„>  ■  0 


n  =  nB  at  a  =  og  ,  r,(ag)  =  0 


*<*> 


and  denoting  the  transformed  metric  tensor  as  g ,  we  have 


(4.2) 


gn  =  guA2  ,  gu  =  Xx  +  yx  +  zx 


g12  -  g12/ex  ,  g12  -  XxXo  +  yxyo  +  ZxZo 


g22  -  g22/e2  ,  g22  -  xj  +  yj  +  Z2 


G3  "  G3/02*2  *  G3  =  &H&22  ~  ^g12^2 


(4.3) 


X*X,Y*Y,Z  =  Z 
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ki  +  V-2  =  +  1«2 


r  =  r/02x2 


Further  noting  that 


hi  ■  <Ixx  -  ^)/i2 


r.  =  r  /0X 

-Cn  'xo 


(4.4) 


r  e 

r  =  (r  ^£_£)/e2 
-nn  ~aa  a 


Using  (4.3)  and  (4.4)  in  Eqs.  (3.14)  -  (3.16),  we  have 


g,„x  -  2g  _x  +  g„x  =  Px  +  Qx  +  XR 

22  xx  12  x°  11  x  o 


(4.5) 


g„„y  -  2g  „y  +  g..y  =  Py  +  Qy  +  YR 

°l2Jxx  12  x°  11  oo  x  a 


(4.6) 


g0„z  -  2g  z  +  g.,z  =  Pz  +  Qz  +  ZR 
22  xx  12  xo  B11  oa  x  o 


(4.7) 


where 


g22 

P  = 

1  X 


811 

Q  =  -rea 


(4.8) 


Thus,  by  choosing  X  and  9  arbitrarily  we  can  redistribute  the  coordinates 
In  the  desired  manner.  An  example  of  this  choice  is  given  in  the  next 
section. 
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5.  An  Analytical  Example  of  Coordinate  Generation 

In  this  section  we  shall  consider  the  problem  of  coordinate  genera¬ 
tion  between  a  prolate  ellipsoid  and  a  sphere  with  coordinate  contraction 
near  the  inner  surface.  This  problem  yields  an  exact  solution  of  the 
equations  (4.5)  -  (4.7). 

Let  n  =  n_  and  n  =  n  be  the  inner  prolate  ellipsoid  and  the  outer 
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sphere  respectively.  The  coordinates  which  vary  on  these  two  surfaces 
are  £  and  £. .We  now  envisage  a  net  of  lines  £  =  const,  and  £  =  const, 
on  these  two  surfaces.  A  curve  on  the  inner  surface  designated  as 
C  =  to  is 

x  =  coshn_cosi; 

Bo  \ 

y  =  sinhn„sinCocos£  )  (5.1) 

z  =  sinhngSini^sinf; 

Similarly,  the  curve  corresponding  to  £ 

x  =  e  cost; 

o 

n»  , 

y  =  e  sinCQcosC 

n« 

z  =  e  sinCQsinC 

Based  on  the  forms  of  the  functions  x,y,z  in  (5.1)  and  (5.2),  we 
assume  the  following  forms  of  x,y,z  for  the  surface  (  = 


C  on  the  outer  surface  is 
o 


(5.2) 
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x  =  f(o)cos£o 

y  =  <|>  (o)sin4ocos£; 

z  =  <t>(o)sin<;  sinE 
o 

The  boundary  conditions  for  f  and  $  are 

f(og)  =  coshnB 

f  Co  )  =  e 

co 

4>  C°B)  =  sinhrig 


<K°„)  =  e 


(5.3) 


(5.4) 


Calculating  the  various  derivatives,  metric  coefficients,  and  all 

other  data  needed  in  the  equations  (4.5)  -  (4.7),  we  get  on  substitution 

an  equation  which  has  sin2CQ  and  cos2^.  Equating  to  zero  the  coefficients 

of  sin2t  and  cos2t  ,  we  obtain 
o  o 


f" 

f' 


11  +  H 
e  <j> 


(5.5) 


41-11  +  11 
•j*'  e  $ 


(5.6) 


where  a  prime  denotes  differentiation  with  respect  to  o.  On  direct  inte¬ 
gration  of  Eqs.  (5.5)  and  (5.6)  under  the  boundary  conditions  (5.4), 
we  get 
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By  taking  a  value  of  K  slightly  greater  than  one  (K  =  1.05  or  1.1),  we 
can  have  sufficient  contraction  of  coordinates  near  the  inner  surface. 

For  the  chosen  problem,  since  the  dependence  on  £  is  simple,  we 
find  that  the  coordinates  between  a  prolate  ellipsoid  and  a  sphere  are 

x  =  [AeBn'0^  +  C]cos^ 

y  =  DeBr^°^  sin£cos£ 

z  =  DeBri^0^sintsin£ 

where  A,  B,  C,  and  D  are  given  in  equation  (5.9). 
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6.  Conclusions 


A  new  method  for  the  generation  of  three-dimensional  coordinates 
between  two  arbitrary  shaped  bodies  has  been  presented.  The  method  is 
based  on  some  simple  differential-geometric  concepts  such  as  the  equations 
of  Gauss  and  the  expressions  for  the  principal  curvatures  of  a  surface. 

The  simplicity  of  the  method  lies  in  solving,  at  one  time,  only  three 
partial  differential  equations  of  the  two-dimensional  type.  This  aspect 
is  bound  to  reduce  the  working  core  requirements  for  a  given  problem  on 
a  computing  machine.  Finally  the  method  allows,  in  a  very  direct  fashion, 
the  possibility  of  coordinate  redistribution  in  the  desired  regions  (cf. 
Eqs.  (4.5)  -  (4.7)). 

An  analytic  solution  of  the  proposed  equations  for  the  case  of  an 
inner  prolate  ellipsoid  and  an  outer  sphere  has  been  presented.  This 
example  shows  that  one  can  generate  coordinates  between  two  analytically 
specified  surfaces  of  simple  forms  by  exact  solutions  of  the  proposed 
equations. 

In  this  paper  the  fundamental  equations  which  form  a  set  of  con¬ 
straints  for  the  generation  of  coordinates  in  the  surface  are 

&2c  =  0 

A2n  «  0  , 

where  A2  is  the  surface  Beltrami  operator  of  the  second  order.  It  must 
be  noted  that  A2  is  neither  a  Laplace  operator  in  the  Cartesian  plane 
(x,y),  nor  in  the  Cartesian  space  (x,y,z).  However,  in  the  case  of  a 
Cartesian  plane  (x,y),  when  there  is  no  dependence  on  z,  A2  reduces  to 
the  Laplace  operator  V2. 
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Introduction 

The  problem  of  generating  a  computational  grid  about  an  aircraft 
configuration  has  recently  been  addressed  by  Ericksson  [4]  and  Lee  [7]. 

In  this  report  the  basic  steps  in  a  grid  generation  scheme  will  be 
described.  No  particular  aircraft  is  considered,  but  the  procedure  is 
intended  to  be  versatile  and  able  to  handle  many  types  of  configurations 
with  only  slight  changes  in  the  computational  space.  The  primary  components 
of  an  aircraft  are  generally  the  wings  and  fuselage.  Therefore,  these 
components  determine  the  basic  structure  of  the  computational  grid. 

The  first,  and  one  of  the  most  difficult,  task  encountered  is  the 
description  of  the  surface  coordinates.  In  the  design  of  various  geometric 
objects,  often  paramatric  surface  patches  are  joined  along  specific  curves 
on  the  surface  of  the  object.  These  patches  may  be  familiar  quadratic 
surfaces,  which  will  be  considered  here,  or  surfaces  defined  by  interpolating 
polynomials  or  splines.  Of  course,  the  shape  of  these  surface  patches 
will  depend  on  the  type  of  grid  which  is  to  surround  the  body.  The  surface 
parametrization  induces  a  natural  grid  system  on  the  body  which  may  not 
be  best  for  computational  purposes.  Thus  a  method  for  reparameterizing 
the  surface  patch  will  be  described.  This  leads  to  considerable  control 
over  the  distribution  of  surface  grid  points. 

For  the  general  configuration  being  considered,  a  single  rectangular 
computational  region  will  be  used.  A  better  distribution  of  grid  points 
can  often  be  obtained  by  adding  or  subtracting  rectangular  blocks  from 
the  computational  region  as  noted  by  Lee  et  al .  [7].  However,  the  basic 
topology  of  the  grid  structure  remains  the  same.  Slits  in  the  computational 
region  may  be  used  to  represent  the  horizontal  and  vertical  stabilizers 


and  any  other  fin-like  components.  The  nacelles  may  be  incorporated  by 
removing  rectangular  blocks  from  the  computational  region. 

The  physical  region  about  the  aircraft,  from  a  topological  viewpoint, 
is  obtained  by  folding  the  computational  region  about  the  body.  This  folding 
results  in  mesh  points  with  irregular  neighborhood  structures.  Although 
the  transformation  from  the  computational  to  the  physical  region  is  singular 
at  such  points,  no  difficulty  is  encountered  in  deriving  accurate  difference 
equations  for  approximating  a  partial  differential  equation. 

The  grid  generation  schemes  discussed  here  are  algebraic  with  emphasis 
on  the  control  of  the  grid  lines.  The  values  of  the  grid  point  coordinates 
can  be  used  as  initial  data  in  an  elliptic  system  which  would  smooth  out 
any  discontinuities  in  grid  line  tangents  or  ripples  caused  by  the  body 
or  interpolating  functions.  No  particular  interpolation  procedure  will 
be  stressed.  Rather,  a  general  format  will  be  followed  which  permits 
inclusion  of  polynomial  interpolation,  splines,  and  the  multi-surface 
method  developed  by  Eiseman  [2,3]. 

Coordinate  System  Defined  by  the  Grid 

The  arrangement  of  the  grid  lines  can  be  visualized  by  considering 
a  curvilinear  coordinate  system  about  the  aircraft  body.  For  fluid  flow 
calculations  it  is  desirable  to  have  the  majority  of  the  grid  points  located 
near  the  body  and  in  the  wake  region.  A  curvilinear  coordinate  system  may 
be  generated  by  folding  a  rectangular  region  about  the  body.  For  example, 
a  first  fold  can  be  made  around  the  front  of  the  body  with  edges  coinciding 
with  the  wing  tips  and  proceeding  downstream.  Folds  at  each  wing  tip  then 
result  in  a  region  enclosing  the  aircraft  as  illustrated  in  Figure  1.  The 
lines  from  the  aircraft  to  the  outer  boundary  are  lines  along  which  the 
coordinate  system  is  singular.  A  more  detailed  discussion  of  this  type 
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of  coordinate  system  has  been  given  by  Ericksson  [4].  This  construction 
results  in  an  ellipsoidal-cylindrical  coordinate  system  surrounding  the 
body.  The  elliptical  portion  gives  good  resolution  of  the  area  around 
the  body  while  the  cylindrical  portion  is  used  to  resolve  the  wake  region. 

Surface  Grid  Points 

The  coordinates  of  the  grid  points  on  the  surface  of  the  aircraft 
and  on  the  outer  boundary  must  be  compatible  with  the  curvilinear  coordinate 
system  surrounding  the  aircraft.  This  will,  be  illustrated  by  considering 
a  simple  wing-fuselage  conf iguration.  All  components  are  defined  using 
linear  and  quadratic  equations  and  the  shape  of  the  body  can  be  changed 
by  redefining  only  a  few  basic  parameters.  Two  particular  variations  are 
given  in  Figures  2  and  3.  The  outer  boundary  is  depicted  in  Figure  4. 

Even  though  the  surface  components  are  relatively  simple,  one  must 
be  able  to  select  the  grid  points  according  to  the  requirements  of  the 
problem  being  solved.  For  example,  one  may  wish  to  have  a  uniform  distri¬ 
bution  of  mesh  points  or  to  have  a  concentration  of  mesh  points  in  regions 
of  particular  interest.  The  general  problem  of  determining  grid  points 
on  an  arbitrary  surface  given  by  parametric  equations  will  now  be  addressed. 

Suppose  a  surface  is  given  by  the  parametric  equation 

V  =  V(s,t) 

where  V  =  (x,y,z)  and  0  <  s,t  <  1.  The  choice  of  a  finite  number  of  s 
and  t  values  will  determine  the  grid  points  on  the  surface.  This  choice 
will  depend  on  the  desired  distribution  of  grid  points.  Although  inter¬ 
polation  on  a  table  of  values  for  s  and  t  could  be  used,  the  following 
procedure  is  very  simple  and  effective.  We  first  consider  the  choice  of 
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s  value  for  a  fixed  t.  A  change  of  variables  to  a  new  parameter  £,  with 
s  =  s (5) »  s'  (?)  >  0,  will  be  assumed  so  that  equally  spaced  5  values  can 
be  used  in  the  computational  region.  The  distribution  of  grid  points  along 
the  t  =  constant  curve  can  be  specified  by  the  arc  length  derivative 

*U)  =  |v5(s(d.t)|  =  |Vs(s,t)|s‘(0. 

Thus  s(c)  must  satisfy  the  differential  equation 

S  U)  =  *U)|vs(s,t)|-1  0) 


which  can  be  solved  by  any  standard  numerical  method  once  an  initial  condition 
is  specified.  In  the  numerical  solution,  the  value  £(c)  is  approximated 
by  the  desired  spacing  between  mesh  points.  The  following  two-dimensional 
example  will  clarify  the  general  procedure.  Let  a  quadrant  of  an  ellipse 
be  given  by 

x  =  a  cos  s,  y  =  b  sin  s,  0  >  s  > 


suppose  equally  spaced  points  are  desired.  The  length  of  the  curve  is 
approximately 


L-if 


a  +  b 


Thus  the  length  of  the  arc  between  any  two  grid  points  is  approximately 
K  s  L/(N-1)  where  N  is  the  number  of  grid  points.  Equation  (1)  then  becomes 


s'(0  -  K[a2  sin2  (s(0)  +  b2  cos2(s 5  S  <  , 


5 


with  initial  condition  s(C0)  =  0.  Due  to  the  error  in  the  numerical  solution 
and  the  approximation  of  arc  length,  some  scaling  was  necessary  in  order 
that  s (c-j )  =  jj-  .  Even  with  this  rescaling,  the  uniform  distribution  of 
grid  points  along  the  central  vertical  coordinate  line  in  Figure  4  attests 
to  the  effectiveness  of  this  technique.  This  ellipse  was  used  to  determine 
one  of  the  parameter  distributions  for  the  ellipsoidal  outer  boundary. 

Equal  spaced  values  for  the  other  parameter  were  used  since  this  automatically 
concentrated  grid  points  near  the  wing  tips. 

The  grid  points  on  the  fuselage  determined  the  position  of  one  set 
of  grid  lines  on  the  wings.  The  other  set  was  chosen  to  concentrate  grid 
points  near  the  wing  tips.  The  same  procedure  would  be  used  for  any  other 
fin-like  structure  projecting  from  the  fuselage. 

Flow  Field  Grid 

For  regions  with  smooth  boundaries,  interpolation  has  been  used  exten¬ 
sively  for  grid  generation  in  two  and  three-dimensional  regions.  Given 
a  set  of  points  ,  ?£,  ....  Pn,  and  a  set  of  parameter  values  s-| ,  $2» 

.  .  .  ,  sR,  a  curve  passing  through  the  points  P-,  with  V(s^)  =  P.-  can 
be  represented  by 

n 

V(s)  =  £  *.(s)P. 

i  =  1  1  1 

where  4.  is  a  function  satisfying  Ms.)  =  6..  . 

1  *  J  '  J 

While  Lagrange  polynomials  could  be  used  for  the  4^,  they  frequently  lead 
to  curves  with  extreme  oscillation  and  cubic  splines  are  a  more  commonly 
used  alternative.  When  the  4..  are  cubic  splines,  their  derivatives  at 
$1  and  sn  may  be  given  arbitrarily.  Thus  the  direction  of  the  tangent 
vector  to  the  curve  at  and  Pn  can  be  specified.  Cubic  splines  have 


continuous  second  order  derivatives  which  produces  curves  sufficiently 
smooth  for  most  grid  generation  problems.  In  the  following  applications 
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the  points  ,  i  =  2,  .  .  .  ,  n  -  1 ,  often  do  not  lie  on  the  body  and  are 
used  only  to  control  the  general  direction  of  the  curve.  In  this  case 
the  condition  on  the  above  can  be  replaced  by 


V*j> 


6  ■  •  for  j  =  1  or  n 
■  J 


n 

and  5Z  fc(s)  =  1 ,  «.  .(s)  >  0 
i  =  1  1  1 


Examples  of  functions  satisfying  these  conditions  are  the  B-spline  basis 
functions  developed  by  de  Boor  [1]  and  used  by  Gordon  and  Riesenfeld  [6] 
and  the  basis  functions  used  by  Eiseman  [3]  in  his  multisurface  method. 

The  latter  development  is  more  general  and  could  include  basis  functions 
other  than  piecewise  polynomials.  It  was  observed  that  the  piecewise 
quadratic  curves  of  Eiseman  are  the  same  as  would  be  generated  using 
B-splines. 

'Returning  to  the  general  configuration  in  Figure  1,  all  that  is  presently 

given  is  the  grid  points  on  the  aircraft  surface  and  on  the  ellipsoidal 

outer  boundary.  Next  a  surface  of  grid  points  will  be  constructed  in  the 

wake  region  behind  the  aircraft.  This  will  employ  interpolating  functions. 

A  set  of  intersecting  curves  is  constructed  as  in  Figure  5  to  connect  the 

aircraft  with  the  downstream  boundary.  These  are  interpolation  curves 

which  are  parameterized  using  surface  parameters  s  and  t  so  that 

r.  »  V(s,t.)  and  A.  =  V(s.,t).  Once  the  basis  functions  are  selected  the 
1  *  J  J 

two-dimensional  blended  interpolant  of  Gordon  and  Hall  [5]  completes  the 
generation  of  the  grid  on  the  wake  surface.  Such  a  surface  is  shown  in 
Figure  6. 


The  grid  points  on  two  coordinate  surfaces,  say  r  =  0  and  r  =  1,  have 
been  determined.  Now  the  space  between  the  two  surfaces,  the  aircraft- 
wake  surface  and  the  outer  boundary,  must  be  filled  with  grid  points.  This 
can  also  be  accomplished  by  interpolation.  Two  intermediate  surfaces  are 
used  to  control  the  interpolating  curves.  The  first  consists  of  points 
at  a  fixed  distance  along  normals  to  the  r  =  0  surface.  The  second  surface 
contains  points  on  the  line  segment  between  these  points  on  the  normal 
surface  and  the  corresponding  points  on  the  outer  surface.  Thus  the  interior 
interpolation  grid  points  are  as  depicted  in  Figure  7.  Figure  8  contains 
plots  of  grids  constructed  about  a  wing-fuselage  configuration.  The  distribu¬ 
tion  of  grid  points  is  controlled  by  the  location  of  the  two  intermediate 
surfaces  and  the  number  of  coordinate  surfaces  between  each  pair  of  inter¬ 
mediate  and  boundary  surfaces.  The  first  surface  must  be  relatively  close 
to  the  aircraft  to  prevent  grid  lines  connecting  the  aircraft  to  the  outer 
boundary  from  crossing. 

The  addition  of  the  other  components  to  the  basic  wing-fuselage 
configuration  necessitates  extra  control  surfaces  in  the  interpolation 
process.  Since  some  of  these  control  surfaces  extend  from  the  aircraft 
to  the  outer  boundary,  a  three-dimensional  blended  interpolant  will  be 
needed  (see  [5]).  The  control  surface  which  would  be  constructed  for  a 
stabilizer  is  depicted  in  Figure  9.  The  construction  of  this  surface  would 
essentially  duplicate  the  previous  construction  of  the  surface  in  the  wake 
region.  Note  that  computationally  both  surfaces  of  the  stabilizer  lie 
on  the  same  control  surface.  The  number  of  nodal  points  used  in  the 
interpolation  scheme  for  constructing  the  grid  above  and  below  the 
stabilizers  would  be  the  same.  However  coordinate  values  on  the  upper 
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surface  would  be  used  in  one  case  while  coordinate  values  on  the  lower 
surface  would  be  used  in  the  other  case. 

Coordinate  Singularities 

It  was  noted  previously  that  the  curvilinear  coordinate  system  about 
the  body  has  singular  lines  along  which  the  Jacobian  of  the  transformation 
from  rectangular  to  curvilinear  coordinates  vanishes.  These  lines  were 
indicated  in  Figure  1.  Since  these  points  have  only  four  immediate 
neighboring  grid  points,  rather  than  the  usual  six,  it  is  probably  best 
to  delete  these  points  from  the  computational  process  when  solving  fluid 
flow  problems.  If  this  is  done  the  grid  structure  about  a  singular  point 
appears  as  in  Figure  10.  Only  the  points  on  the  surface  cutting  the 
singular  line  are  shown.  Now  all  first  and  second  order  difference 
expressions  can  be  evaluated  in  the  usual  manner.  This  procedure  is 
certainly  simpler  than  using  series  expansions  to  derive  separate  differ¬ 
ence  equations  at  the  singular  points.  Previous  calculations  by  Mastin, 
Ghosh,  and  Thompson  [8]  have  shown  that  this  technique  produces  accurate 
results  in  potential  and  viscous  flow  problems. 
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Figure  1.  Aircraft  and  flow  field 


Figure  2.  Wing-fuselage  configuration  defined  by 
elementary  surfaces. 


Figure  4.  Outer  boundary  of  flow  field. 
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Figure  5.  Interpolation  curves  for  wake  surface 


Figure  7. 


Interpolation  points  between  aircraft 
surface  and  outer  boundary. 


Grid  on  interpolati 
for  stabilizer. 


